%% Prepare Workspace and Load Data
delete(gcp('nocreate')) %Close the Parallel Pools
close all; clear all; clc; %Clear workspace
load Sall_SallCr_151w_23.mat; %Load Endmember Matrix
load ImageNames.mat; %Load List of Image Names

%%
for k = 15 %loop through all 17 images
    parpool(10) %Initalize Parallel Pool
    name = ImageNames(k); %Select kth image name
    c = clock; %Note the start time
    cd 'F:\FinalMaps\2 Cropped Raw Images'
    % Load/Read kth Image
    Hdr = readlines(append(name,"_10000_int_179.hdr")); %read kth image header
    NS = sscanf(Hdr(4), 'samples = %f'); %extract number of samples from image header
    NL = sscanf(Hdr(5), 'lines = %f'); %extract number of lines from image header
    NB = sscanf(Hdr(6), 'bands = %f'); %extract number of bands from image header
    img = multibandread(append(name,"_10000_int_179.img"),[NL NS NB],'int16',0,'bil','ieee-le'); %load kth image
    img(:,:,[70:93,176:179]) = []; % crop image to wavelength region of interest
    
    % Run through a Routine to NaN Dark Pixels
    img = ((sum(img,3) > 180000)&(sum(img,3) > 0)).*img; %zero out any dark pixels with a total summed reflectance below 240000 (determined empirically)
    img(img==0) = NaN; %replace these zeros with NaNs
    sizeimg = size(img); %capture image size
    rows = sizeimg(1); %capture number of rows of image
    columns = sizeimg(2) ;%capture number of columns of image
    
    %Initialize Matrices for Output of Misfit Map, Final Abundance
    %(Model Parameter) Maps, and the Model Prediction
    Fmstore = NaN(sizeimg(1),sizeimg(2),1); %Initialize store matrix for map of model misfit at each pixel across the image
    abundances = NaN(sizeimg(1),sizeimg(2),width(Sall)); %Initialize a store matrix for the end member abundance values (i.e., renomalized model parameters) at each pixel across the image
    modelprediction = NaN(size(img)); %Initialize the store matrix for the modelled spectrum at each pixel across the image
    
    %Unmix Current End Members and Store Associated Products
    parfor i = 1:rows %loop through all rows of image
        for j = 1:columns %loop through all columns of image
                d = squeeze(img(i,j,:)); % select spectrum for individual pixel
            if  isnan(d(1)) %if this spectrum is a NaN
                Fmstore(i,j,1) = NaN; %store a misfit value of a NaN
                abundances(i,j,:)= NaN; %store abundance values of a NaN
                modelprediction(i,j,:) = NaN; %store model prediction values of a NaN
            else %otherwise
                % Remove Linear Continuum From Pixel Spectrum
                Gcr = [ones(length(w),1) w]; %Create G matrix for a line through the spectrum
                m = lsqlin(Gcr,d); %Calculate slope & y int of best fit line through the spectrum
                gmc = m(1) + w*(m(2)); %Calculate line model prediction vector
                dCr = d - (m(1) + w*(m(2)));  %Remove line from spectrum data
    
                %Unmix End Members from CR Pixel Spectrum
                m = lsqnonneg(SallCr,dCr); %perform non negative unmixing
                gmn = 0; %intilize modelled spectrum
                for p = 1:length(m) %loop through each model parameter, multiply it by it's end member spectrum, and sum across all end members
                    gmni = m(p)*SallCr(:,p);
                    gmn = gmn + gmni; %arrive at final modelled spectrum
                end
                g = gmn + gmc; %Add back in the continuum removed portion of the spectrum 
                Fmstore(i,j,1) =(1/2)*(d-g)'*(d-g); %calculated and store the L2 misfit for the pixel
                abundances(i,j,:) = squeeze(m/sum(m)).*100; %store spectral contributon abundances per end member
                modelprediction(i,j,:) = g; % store model prediction
            end
        end
    end        
    % Create Final Map by Selecting Maximum Renormalized Abundance 
    final_map = zeros(NL,NS);
    for i=1:NL
        for j = 1:NS   
            if sum(isnan(a100(i,j,:))) == 13
                final_map(i,j) = NaN;      
            else
                [val,final_map(i,j)] = max(squeeze(a100(i,j,:)));     
            end
        end
    end

    % Isolate Rock Type Spectral Contributions / Abundances
    a1 = abundances(:,:,1); % gyps 1 gypsym evaporite
    a2 = abundances(:,:,2); % gyps 2 gypsum cemented ss
    a3 = abundances(:,:,3); % gyps 3 clay bearing gypsum evaporite
    a4 = abundances(:,:,4); % dickite ss
    a5 = abundances(:,:,5); % alluvium
    a6 = abundances(:,:,6); % mixed carbonate & ss
    a7 = abundances(:,:,7); % kaolinte ss
    a8 = abundances(:,:,8); % dolomite
    a9 = abundances(:,:,9); % reef limestone out of place
    a10 = abundances(:,:,10); % reef limestone in place
    a11 = abundances(:,:,11); % flat limestone
    a12 = abundances(:,:,12); % dolomite veg 1
    a14 = abundances(:,:,14); % dolomite veg 2

    % Apply Emprical Cutoffs to Each Rock Type End Member
    a1(a1<7) = NaN; % gyps 1 gypsym evaporite
    a2(a2<12) = NaN; % gyps 2 gypsum cemented ss
    a3(a3<10) = NaN; % gyps 3 clay bearing gypsum evaporite
    a4(a4<4) = NaN; % dickite sandstone
    a5(a5<14) = NaN; % alluvium
    a6(a6<7) = NaN; % mixed carbonate & kaolinite sandstone
    a7(a7<15) = NaN; % kaolinte ss
    a8(a8<8) = NaN; % dolomite
    a9(a9<0) = NaN; % reef limestone out of place
    a10(a10<0) = NaN; % reef limestone in place
    a11(a11<4) = NaN; % spectrally flat limestone
    a12(a12<20) = NaN; % dolomite + veg 1
    a14(a14<20) = NaN; % dolomite + veg 2
    
    % Renormalize 13 Rock Type End Member Spectral Contributions to Sum to 100
    abund = cat(3,a1,a2,a3,a5,a4,a7,a6,a12,a14,a8,a9,a10,a11); %Re order to be consistent with manuscript ordering
    abund(isnan(abund)) = 0; %set NaNs to 0
    at = zeros(NL,NS); %initialize summed abundances store matrix
    a100 = zeros(size(abund)); %initialize renormalized abundances store matrix
    for i=1:NL %loop through all image lines
        for j = 1:NS %loop through all image samples
            at(i,j) = sum(squeeze(abund(i,j,:))); %calculate summed abundance for pixel
            for p = 1:13 %for each end member
                a100(i,j,p) = abund(i,j,p)./at(i,j).*100; %calculate renormalized abundance
            end
        end
    end
    a100(a100==0) = NaN; %set 0s to NaN
  
    % Save ENVI Images & Headers
    offset = 0; % Set offset to 0 for ENVI Header
    datatype = 2;  % Set data type for 16 bit signed integer
    byteorder = 0; % Set byte order to Host (intel)
        
    % Model Prediction ENVI Image & Header Save
    cd 'F:\FinalMaps\3 Model Prediction Images' 
    multibandwrite(single(modelprediction),append(name,"_10000_int_151_model_pred.img"),'bil','machfmt','ieee-le','OFFSET',0,'PRECISION','int16')
    [NL,NS,NB] = size(modelprediction);
    fid = fopen(append(name,"_10000_int_151_model_pred.hdr"),'wt');
    text = 'ENVI'; fprintf(fid,'%s\n',text);
    text = sprintf('description = {'); fprintf(fid,'%s\n',text);
    text = sprintf('File Imported into ENVI.}'); fprintf(fid,'%s\n',text);
    text = sprintf('samples = %d',NS); fprintf(fid,'%s\n',text);
    text = sprintf('lines = %d',NL); fprintf(fid,'%s\n',text);
    text = sprintf('bands = %d',NB); fprintf(fid,'%s\n',text);
    text = sprintf('header offset = %d',offset); fprintf(fid,'%s\n',text);
    text = sprintf('file type = ENVI'); fprintf(fid,'%s\n',text);
    text = sprintf('data type = %d',datatype); fprintf(fid,'%s\n',text);
    text = sprintf('interleave = bil'); fprintf(fid,'%s\n',text);
    text = sprintf('sensor type = Unknown'); fprintf(fid,'%s\n',text);
    text = sprintf('byte order = %d',byteorder); fprintf(fid,'%s\n',text);
    text = sprintf(Hdr(13)); fprintf(fid,'%s\n',text);
    text = sprintf('wavelength units = Nanometers'); fprintf(fid,'%s\n',text);
    text = sprintf('wavelength = {1433.549576, 1438.559576, 1443.569576, 1448.579576, 1453.579576, 1458.589576, 1463.599576, 1468.609576, 1473.619576, 1478.629576, 1483.639576, 1488.639576, 1493.649576, 1498.659576, 1503.669576, 1508.679576, 1513.689576, 1518.699576, 1523.709576, 1528.709576, 1533.719576, 1538.729576, 1543.739576, 1548.749576, 1553.759576, 1558.769576, 1563.769576, 1568.779576, 1573.789576, 1578.799576, 1583.809576, 1588.819576, 1593.829576, 1598.839576, 1603.839576, 1608.849576, 1613.859576, 1618.869576, 1623.879576, 1628.889576, 1633.899576, 1638.899576, 1643.909576, 1648.919576, 1653.929576, 1658.939576, 1663.949576, 1668.959576, 1673.969576, 1678.969576, 1683.979576, 1688.989576, 1693.999576, 1699.009576, 1704.019576, 1709.029576, 1714.029576, 1719.039576, 1724.049576, 1729.059576, 1734.069576, 1739.079576, 1744.089576, 1749.099576, 1754.099576, 1759.109576, 1764.119576, 1769.129576, 1774.139576, 2074.659576, 2079.669576, 2084.679576, 2089.689576, 2094.689576, 2099.699576, 2104.709576, 2109.719576, 2114.729576, 2119.739576, 2124.749576, 2129.749576, 2134.759576, 2139.769576, 2144.779576, 2149.789576, 2154.799576, 2159.809576, 2164.819576, 2169.819576, 2174.829576, 2179.839576, 2184.849576, 2189.859576, 2194.869576, 2199.879576, 2204.879576, 2209.889576, 2214.899576, 2219.909576, 2224.919576, 2229.929576, 2234.939576, 2239.949576, 2244.949576, 2249.959576, 2254.969576, 2259.979576, 2264.989576, 2269.999576, 2275.009576, 2280.009576, 2285.019576, 2290.029576, 2295.039576, 2300.049576, 2305.059576, 2310.069576, 2315.079576, 2320.079576, 2325.089576, 2330.099576, 2335.109576, 2340.119576, 2345.129576, 2350.139576, 2355.139576, 2360.149576, 2365.159576, 2370.169576, 2375.179576, 2380.189576, 2385.199576, 2390.209576, 2395.209576, 2400.219576, 2405.229576, 2410.239576, 2415.249576, 2420.259576, 2425.269576, 2430.269576, 2435.279576, 2440.289576, 2445.299576, 2450.309576, 2455.319576, 2460.329576, 2465.339576, 2470.339576, 2475.349576, 2480.359576}'); fprintf(fid,'%s\n',text);
    text = sprintf('data ignore value = NaN');
    fprintf(fid,'%s\n',text);
    fclose(fid);
    
    %  Model Misfit ENVI Image & Header Save
    cd 'F:\FinalMaps\5 Model Misfit Images'
    multibandwrite(single(Fmstore),append(name,"_10000_int_151_model_misfit.img"),'bil','machfmt','ieee-le','OFFSET',0,'PRECISION','float32')
    [NL,NS,NB] = size(Fmstore);
    fid = fopen(append(name,"_10000_int_151_model_misfit.hdr"),'wt');
    text = 'ENVI'; fprintf(fid,'%s\n',text);
    text = sprintf('description = {'); fprintf(fid,'%s\n',text);
    text = sprintf('File Imported into ENVI.}'); fprintf(fid,'%s\n',text);
    text = sprintf('samples = %d',NS); fprintf(fid,'%s\n',text);
    text = sprintf('lines = %d',NL); fprintf(fid,'%s\n',text);
    text = sprintf('bands = %d',NB); fprintf(fid,'%s\n',text);
    text = sprintf('header offset = %d',offset);  fprintf(fid,'%s\n',text); 
    text = sprintf('file type = ENVI');  fprintf(fid,'%s\n',text);
    text = sprintf('data type = %d',4);  fprintf(fid,'%s\n',text);
    text = sprintf('interleave = bil');  fprintf(fid,'%s\n',text);
    text = sprintf('sensor type = Unknown');  fprintf(fid,'%s\n',text);
    text = sprintf('byte order = %d',byteorder);  fprintf(fid,'%s\n',text);
    text = sprintf(Hdr(13));  fprintf(fid,'%s\n',text);
    text = sprintf('data ignore value = NaN'); fprintf(fid,'%s\n',text);
    fclose(fid);
    
    % Save Renormalized Spectral Contribution Maps, ENVI Image & Header
    cd 'F:\FinalMaps\6 Abundance Renormalized Images'
    multibandwrite(single(a100),append(name,"_10000_int_151_abundances_renormalized.img"),'bil','machfmt','ieee-le','OFFSET',0,'PRECISION','int16')
    [NL,NS,NB] = size(a100);
    offset = 0;
    datatype = 2;
    byteorder = 0;
    fid = fopen(append(name,"_10000_int_151_abundances_renormalized.hdr"),'wt');
    text = 'ENVI'; fprintf(fid,'%s\n',text);
    text = sprintf('description = {'); fprintf(fid,'%s\n',text);
    text = sprintf('File Imported into ENVI.}'); fprintf(fid,'%s\n',text);
    text = sprintf('samples = %d',NS); fprintf(fid,'%s\n',text);
    text = sprintf('lines = %d',NL); fprintf(fid,'%s\n',text);
    text = sprintf('bands = %d',NB); fprintf(fid,'%s\n',text);
    text = sprintf('header offset = %d',offset); fprintf(fid,'%s\n',text);  
    text = sprintf('file type = ENVI'); fprintf(fid,'%s\n',text); 
    text = sprintf('data type = %d',datatype); fprintf(fid,'%s\n',text);
    text = sprintf('interleave = bil'); fprintf(fid,'%s\n',text);
    text = sprintf('sensor type = Unknown'); fprintf(fid,'%s\n',text);
    text = sprintf('byte order = %d',byteorder); fprintf(fid,'%s\n',text);
    text = sprintf(Hdr(13)); fprintf(fid,'%s\n',text);
    text = sprintf('data ignore value = NaN');
    fprintf(fid,'%s\n',text);
    fclose(fid);
    
    % Save Final Composite Lithology-Hyperspectral Map, ENVI Image & Header
    cd 'F:\FinalMaps\7 Final Map Images'
    multibandwrite(single(final_map),append(name,"_10000_int_151_final_map.img"),'bil','machfmt','ieee-le','OFFSET',0,'PRECISION','int16')
    [NL,NS,NB] = size(final_map);
    fid = fopen(append(name,"_10000_int_151_final_map.hdr"),'wt');
    text = 'ENVI'; fprintf(fid,'%s\n',text);
    text = sprintf('description = {'); fprintf(fid,'%s\n',text);
    text = sprintf('File Imported into ENVI.}'); fprintf(fid,'%s\n',text);
    text = sprintf('samples = %d',NS); fprintf(fid,'%s\n',text);
    text = sprintf('lines = %d',NL); fprintf(fid,'%s\n',text);
    text = sprintf('bands = %d',NB); fprintf(fid,'%s\n',text);
    text = sprintf('header offset = %d',offset); fprintf(fid,'%s\n',text);  
    text = sprintf('file type = ENVI'); fprintf(fid,'%s\n',text); 
    text = sprintf('data type = %d',datatype); fprintf(fid,'%s\n',text);
    text = sprintf('interleave = bil'); fprintf(fid,'%s\n',text);
    text = sprintf('sensor type = Unknown'); fprintf(fid,'%s\n',text);
    text = sprintf('byte order = %d',byteorder); fprintf(fid,'%s\n',text);
    text = sprintf(Hdr(13)); fprintf(fid,'%s\n',text);
    text = sprintf('data ignore value = NaN');
    fprintf(fid,'%s\n',text);
    fclose(fid);
    
    %Clock Time for Display, Close ParPool & Clear workspace
    c2 = clock; %save the end clock time for kth image
    display(append('Image ',num2str(k),' Finished Elapsed time: ',num2str(floor((c2(4)+ (c2(5)./60)) - (c(4)+ (c(5)./60)))),... %calculate and display elapsed time
        ' hours and ',num2str(((c2(4)+ (c2(5)./60)) - (c(4)+ (c(5)./60)- floor((c2(4)+ (c2(5)./60)) - (c(4)+ (c(5)./60)))))*60), ' mins'))
    clearvars -except ImageNames Sall SallCr w w151  %clear large variables
    cd 'F:\FinalMaps\1 Code'
end